{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.008106,
     "end_time": "2021-04-23T10:41:33.911944",
     "exception": false,
     "start_time": "2021-04-23T10:41:33.903838",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Double/Debiased ML for Partially Linear IV Model\n",
    "\n",
    "References: \n",
    "\n",
    "https://arxiv.org/abs/1608.00060\n",
    "\n",
    "\n",
    "https://www.amazon.com/Business-Data-Science-Combining-Accelerate/dp/1260452778\n",
    "\n",
    "The code is based on the book.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.006963,
     "end_time": "2021-04-23T10:41:33.926085",
     "exception": false,
     "start_time": "2021-04-23T10:41:33.919122",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "\n",
    "# Partially Linear IV Model\n",
    "\n",
    "We consider the partially linear structural equation model:\n",
    "\\begin{eqnarray}\n",
    " &  Y - D\\theta_0 = g_0(X) + \\zeta,  & E[\\zeta \\mid Z,X]= 0,\\\\\n",
    "  & Z = m_0(X) +  V,   &  E[V \\mid X] = 0. \n",
    "\\end{eqnarray}\n",
    "\n",
    "\n",
    "Note that this model is not a regression model unless $Z=D$.  The model  is a canonical\n",
    "model in causal inference, going back to P. Wright's work on IV methods for estimaing demand/supply equations, with the modern difference being that $g_0$ and $m_0$ are nonlinear, potentially complicated functions of high-dimensional $X$.  \n",
    "\n",
    "\n",
    "The idea of this model is that there is a structural or causal relation between $Y$ and $D$, captured by $\\theta_0$, and $g_0(X) + \\zeta$ is the stochastic error, partly explained by covariates $X$.  $V$ and $\\zeta$ are stochastic errors that are not explained by $X$. Since $Y$ and $D$ are jointly determined, we need an external factor, commonly referred to as an instrument, $Z$ to create exogenous variation\n",
    "in $D$.   Note that $Z$ should affect $D$.  The $X$ here serve again as confounding factors, so we can think of variation in $Z$ as being exogenous only conditional on $X$. \n",
    "\n",
    "\n",
    "The causal DAG this model corresponds to is given by:\n",
    "$$\n",
    "Z \\to D,  X \\to (Y, Z, D),  L \\to (Y,D),\n",
    "$$\n",
    "where $L$ is the latent confounder affecting both $Y$ and $D$, but not $Z$.\n",
    "\n",
    "\n",
    "\n",
    "---\n",
    "\n",
    "# Example\n",
    "\n",
    "A simple contextual example is from biostatistics, where $Y$ is a health outcome and $D$ is indicator of smoking.  Thus, $\\theta_0$ is captures the effect of smoking on health.  Health outcome $Y$ and smoking behavior $D$ are treated as being jointly determined.  $X$ represents patient characteristics, and $Z$ could be a doctor's advice not to smoke (or another behavioral treatment) that may affect the outcome $Y$ only through shifting the behavior $D$, conditional on characteristics $X$.   \n",
    "\n",
    "----\n",
    "\n",
    "\n",
    "\n",
    "# PLIVM in Residualized Form\n",
    "\n",
    "The PLIV model above can be rewritten in the following residualized form:\n",
    "$$\n",
    "  \\tilde Y = \\tilde D \\theta_0 + \\zeta,   \\quad  E[\\zeta \\mid V,X]= 0,\n",
    "$$\n",
    "where\n",
    "$$\n",
    " \\tilde Y = (Y- \\ell_0(X)),  \\quad \\ell_0(X) = E[Y \\mid X] \\\\\n",
    "   \\tilde D = (D - r_0(X)), \\quad r_0(X) = E[D \\mid X] \\\\\n",
    "   \\tilde Z = (Z- m_0(X)), \\quad m_0(X) = E[Z \\mid X].\n",
    "$$\n",
    "   The tilded variables above represent original variables after taking out or \"partialling out\"\n",
    "  the effect of $X$.  Note that $\\theta_0$ is identified from this equation if $V$ \n",
    "  and $U$ have non-zero correlation, which automatically means that $U$ and $V$\n",
    "  must have non-zero variation.\n",
    "\n",
    "  \n",
    "\n",
    "-----\n",
    "\n",
    "# DML for PLIV Model\n",
    "\n",
    "Given identification, DML  proceeds as follows\n",
    "\n",
    "Compute the estimates $\\hat \\ell_0$, $\\hat r_0$, and $\\hat m_0$ , which amounts\n",
    "to solving the three problems of predicting $Y$, $D$, and $Z$ using\n",
    "$X$, using any generic  ML method, giving us estimated residuals \n",
    "$$\n",
    "\\tilde Y = Y - \\hat \\ell_0(X), \\\\ \\tilde D= D - \\hat r_0(X), \\\\ \\tilde Z = Z- \\hat m_0(X).\n",
    "$$ \n",
    "The estimates should be of a cross-validated form, as detailed in the algorithm below. \n",
    "\n",
    "Estimate $\\theta_0$ by the the intstrumental\n",
    "variable regression of $\\tilde Y$ on $\\tilde D$ using $\\tilde Z$ as an instrument.\n",
    "Use the conventional inference for the IV regression estimator, ignoring\n",
    "the estimation error in these residuals. \n",
    "\n",
    "The reason we work with this residualized form is that it eliminates the bias\n",
    "arising when solving the prediction problem in stage 1. The role of cross-validation\n",
    "is to avoid another source of bias due to potential overfitting.\n",
    "\n",
    "The estimator is adaptive,\n",
    "in the sense that the first stage estimation errors do not affect the second \n",
    "stage errors.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "_kg_hide-output": true,
    "papermill": {
     "duration": 34.670095,
     "end_time": "2021-04-23T10:42:08.603197",
     "exception": false,
     "start_time": "2021-04-23T10:41:33.933102",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "install.packages(\"hdm\")\n",
    "install.packages(\"AER\")\n",
    "install.packages(\"randomForest\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 1.846109,
     "end_time": "2021-04-23T10:42:10.458406",
     "exception": false,
     "start_time": "2021-04-23T10:42:08.612297",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "\n",
    "library(AER)  #applied econometrics library\n",
    "library(randomForest)  #random Forest library\n",
    "library(hdm) #high-dimensional econometrics library\n",
    "library(glmnet) #glm net\n",
    "\n",
    "\n",
    "# DML for PLIVM\n",
    "\n",
    "DML2.for.PLIVM <- function(x, d, z, y, dreg, yreg, zreg, nfold=2) {\n",
    "  # this implements DML2 algorithm, where there moments are estimated via DML, before constructing\n",
    "  # the pooled estimate of theta randomly split data into folds\n",
    "  nobs <- nrow(x)\n",
    "  foldid <- rep.int(1:nfold,times = ceiling(nobs/nfold))[sample.int(nobs)]\n",
    "  I <- split(1:nobs, foldid)\n",
    "  # create residualized objects to fill\n",
    "  ytil <- dtil <- ztil<- rep(NA, nobs)\n",
    "  # obtain cross-fitted residuals\n",
    "  cat(\"fold: \")\n",
    "  for(b in 1:length(I)){\n",
    "    dfit <- dreg(x[-I[[b]],], d[-I[[b]]])  #take a fold out\n",
    "    zfit <- zreg(x[-I[[b]],], z[-I[[b]]])  #take a fold out\n",
    "    yfit <- yreg(x[-I[[b]],], y[-I[[b]]])  # take a folot out\n",
    "    dhat <- predict(dfit, x[I[[b]],], type=\"response\")  #predict the fold out\n",
    "    zhat <- predict(zfit, x[I[[b]],], type=\"response\")  #predict the fold out\n",
    "    yhat <- predict(yfit, x[I[[b]],], type=\"response\")  #predict the fold out\n",
    "    dtil[I[[b]]] <- (d[I[[b]]] - dhat) #record residual\n",
    "    ztil[I[[b]]] <- (z[I[[b]]] - zhat) #record residual\n",
    "    ytil[I[[b]]] <- (y[I[[b]]] - yhat) #record residial\n",
    "    cat(b,\" \")\n",
    "  }\n",
    "  ivfit= tsls(y=ytil,d=dtil, x=NULL, z=ztil, intercept=FALSE)\n",
    "  coef.est <-  ivfit$coef          #extract coefficient \n",
    "  se <-  ivfit$se                  #record standard error\n",
    "  cat(sprintf(\"\\ncoef (se) = %g (%g)\\n\", coef.est , se))\n",
    "  return( list(coef.est =coef.est , se=se, dtil=dtil, ytil=ytil, ztil=ztil) )\n",
    "}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.011698,
     "end_time": "2021-04-23T10:42:10.482689",
     "exception": false,
     "start_time": "2021-04-23T10:42:10.470991",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "-----\n",
    "\n",
    "# Emprical Example:  Acemoglu, Jonsohn, Robinson (AER).\n",
    "\n",
    "\n",
    "* Y is log GDP;\n",
    "* D is a measure of Protection from Expropriation, a proxy for \n",
    "quality of insitutions;\n",
    "* Z is the log of Settler's mortality;\n",
    "* W are geographical variables (latitude, latitude squared, continent dummies as well as interactions)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 5.743953,
     "end_time": "2021-04-23T10:42:16.240693",
     "exception": false,
     "start_time": "2021-04-23T10:42:10.496740",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "\n",
    "data(AJR); \n",
    "\n",
    "y = AJR$GDP; \n",
    "d = AJR$Exprop; \n",
    "z = AJR$logMort\n",
    "xraw= model.matrix(~ Latitude+ Africa+Asia + Namer + Samer, data=AJR)\n",
    "x = model.matrix(~ -1 + (Latitude + Latitude2 + Africa + \n",
    "                           Asia + Namer + Samer)^2, data=AJR)\n",
    "dim(x)\n",
    "\n",
    "# DML with Random Forest\n",
    "cat(sprintf(\"\\n DML with Random Forest \\n\"))\n",
    "\n",
    "dreg <- function(x,d){ randomForest(x, d) }  #ML method=Forest\n",
    "yreg <- function(x,y){ randomForest(x, y) }  #ML method=Forest\n",
    "zreg<-  function(x,z){ randomForest(x, z)}  #ML method=Forest \n",
    "  \n",
    "set.seed(1)\n",
    "DML2.RF = DML2.for.PLIVM(xraw, d, z, y, dreg, yreg, zreg, nfold=20)\n",
    "\n",
    "# DML with PostLasso\n",
    "cat(sprintf(\"\\n DML with Post-Lasso \\n\"))\n",
    "\n",
    "dreg <- function(x,d){ rlasso(x, d) }  #ML method=lasso\n",
    "yreg <- function(x,y){ rlasso(x, y) }  #ML method=lasso\n",
    "zreg<-  function(x,z){ rlasso(x, z)}  #ML method=lasso \n",
    "\n",
    "set.seed(1)\n",
    "DML2.lasso = DML2.for.PLIVM(x, d, z, y, dreg, yreg, zreg, nfold=20)\n",
    "\n",
    "\n",
    "# Compare Forest vs Lasso\n",
    "\n",
    "comp.tab= matrix(NA, 3, 2)\n",
    "comp.tab[1,] = c( sqrt(mean((DML2.RF$ytil)^2)),  sqrt(mean((DML2.lasso$ytil)^2)) )\n",
    "comp.tab[2,] = c( sqrt(mean((DML2.RF$dtil)^2)),  sqrt(mean((DML2.lasso$dtil)^2)) )\n",
    "comp.tab[3,] = c( sqrt(mean((DML2.RF$ztil)^2)),  sqrt(mean((DML2.lasso$ztil)^2)) )\n",
    "rownames(comp.tab) = c(\"RMSE for Y:\", \"RMSE for D:\", \"RMSE for Z:\")\n",
    "colnames(comp.tab) = c(\"RF\", \"LASSO\")\n",
    "print(comp.tab, digits=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.013703,
     "end_time": "2021-04-23T10:42:16.269251",
     "exception": false,
     "start_time": "2021-04-23T10:42:16.255548",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Examine if we have weak instruments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 16.107321,
     "end_time": "2021-04-23T10:42:32.390552",
     "exception": false,
     "start_time": "2021-04-23T10:42:16.283231",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "install.packages(\"lfe\")\n",
    "library(lfe)\n",
    "summary(felm(DML2.lasso$dtil~DML2.lasso$ztil), robust=T)\n",
    "summary(felm(DML2.RF$dtil~DML2.RF$ztil), robust=T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.015919,
     "end_time": "2021-04-23T10:42:32.423865",
     "exception": false,
     "start_time": "2021-04-23T10:42:32.407946",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "#  We do have weak instruments, because t-stats in regression $\\tilde D \\sim \\tilde Z$ are less than 4 in absolute value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.015943,
     "end_time": "2021-04-23T10:42:32.455719",
     "exception": false,
     "start_time": "2021-04-23T10:42:32.439776",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "So let's carry out DML inference combined with Anderson-Rubin Idea"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.036055,
     "end_time": "2021-04-23T10:42:32.507790",
     "exception": false,
     "start_time": "2021-04-23T10:42:32.471735",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# DML-AR (DML with Anderson-Rubin) \n",
    "\n",
    "DML.AR.PLIV<- function(rY, rD, rZ, grid, alpha=.05){\n",
    "    n = length(rY)\n",
    "    Cstat = rep(0, length(grid))\n",
    "    for (i in 1:length(grid)) {\n",
    "    Cstat[i] <-  n* (mean( (rY - grid[i]*rD)*rZ)  )^2/var ( (rY - grid[i]*rD) * rZ )\n",
    "    };\n",
    "    LB<- min(grid[ Cstat < qchisq(1-alpha,1)]);\n",
    "    UB <- max(grid[ Cstat < qchisq(1-alpha,1)]); \n",
    "    plot(range(grid),range(c( Cstat)) , type=\"n\",xlab=\"Effect of institutions\", ylab=\"Statistic\", main=\" \");  \n",
    "    lines(grid, Cstat,   lty = 1, col = 1);       \n",
    "    abline(h=qchisq(1-alpha,1), lty = 3, col = 4);\n",
    "    abline(v=LB, lty = 3, col = 2);\n",
    "    abline(v=UB, lty = 3, col = 2);\n",
    "    return(list(UB=UB, LB=LB))\n",
    "    }\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.478528,
     "end_time": "2021-04-23T10:42:33.002469",
     "exception": false,
     "start_time": "2021-04-23T10:42:32.523941",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "DML.AR.PLIV(rY = DML2.lasso$ytil, rD= DML2.lasso$dtil, rZ= DML2.lasso$ztil,\n",
    "           grid = seq(-2, 2, by =.01))\n",
    "\n",
    "\n",
    "DML.AR.PLIV(rY = DML2.RF$ytil, rD= DML2.RF$dtil, rZ= DML2.RF$ztil,\n",
    "           grid = seq(-2, 2, by =.01))\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.3"
  },
  "papermill": {
   "duration": 62.79319,
   "end_time": "2021-04-23T10:42:34.038582",
   "environment_variables": {},
   "exception": null,
   "input_path": "__notebook__.ipynb",
   "output_path": "__notebook__.ipynb",
   "parameters": {},
   "start_time": "2021-04-23T10:41:31.245392",
   "version": "2.1.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
